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Abstract 


The  objective  of  this  document  is  to  investigate  the  feasibility  of  using  the  LS-DYNA  code  to 
perform  a  finite  element  analysis  to  study  the  structural  damage  of  a  typical  compartment-type 
vessel  caused  by  an  internal  blast.  The  work  is  divided  into  two  investigations:  one  considers  the 
generation  and  propagation  of  the  blast  and  the  other  investigation  looks  at  the  material-damage 
and  failure  models  that  could  be  used  to  predict  the  damage  on  the  compartment  structure. 

The  first  investigation  presents  results  and  lessons  learned  from  four  studies.  These  studies 
simulated  the  blast  propagation  using  the  arbitrary  lagrangian  eulerian  (ALE)  approach  in  LS- 
DYNA.  This  investigation  included  parametric  studies  as  well.  Regardless  of  the  work  done,  it 
was  not  possible  to  generate  a  finite  element  model,  fine  enough  (and  yet  still  manageable)  to 
capture  maximum  incident  pressure  and  impulse  with  a  domain  large  enough  to  include  the  whole 
compartment.  Alternative  methods,  such  as  the  raytracer  approach,  were  not  part  of  this  work  but 
should  be  addressed  in  the  future. 

The  second  investigation  discusses  different  issues  regarding  the  use  of  damage  and  failure 
models  in  LS-DYNA.  It  is  recommended  to  develop  an  ad  hoc  method  to  predict  the  damage  and 
failure  in  large  scale  model.  This  method  should  include  a  material- damage  model,  a  failure 
criterion,  as  well  as  a  study  on  the  mesh.  Experimental  data  should  also  be  used  to  validate  the 
model. 
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Resume 


L’objectif  de  ce  document  est  d’explorer  la  faisabilite  d’utiliser  le  programme  LS-DYNA  pour 
realiser  une  analyse  par  elements  finis  permettant  d’etudier  l’endommagement  structurel  d’un 
compartiment  typique  d’un  bateau  cause  par  un  effet  de  souffle  interne.  Le  travail  est  divise  en 
deux  investigations  principals  :  l’une  pour  generer  et  propager  1’ effet  de  souffle  et  1’ autre 
investigation  pour  explorer  les  modeles  d’endommagement  et  de  rupture  qui  pourraient  etre 
utilises  pour  predire  l’endommagement  cause  par  une  explosion  interne  sur  la  structure  du 
compartiment. 

La  premiere  partie  du  travail  presente  les  resultats  et  les  lecons  apprises  pour  quatre  etudes.  Ces 
etudes  ont  simule  la  propagation  de  l’explosion  en  utilisant  l’approche  arbitraire  lagrangienne 
eulerienne  (ALE)  dans  LS-DYNA.  Cette  investigation  incluait  egalement  plusieurs  etudes 
parametriques.  Malgre  tout  le  travail  realise,  il  n’a  pas  ete  possible  de  generer  un  modele  elements 
finis,  suffisamment  raffine  (tout  en  restant  gerable)  pour  predire  adequatement  la  pression 
maximale  incidente  et  1’ impulsion  dans  un  domaine  suffisamment  large  pour  inclure  le 
compartiment  au  complet.  Des  methodes  alternatives,  telle  que  l’approche  ‘raytracer’,  n’ont  pas 
fait  l’objet  de  ce  travail,  mais  devraient  etre  etudiees  dans  le  futur. 

La  deuxieme  partie  du  travail  discute  des  differentes  voies  concernant  l’utilisation  des  modeles 
d’endommagement  et  de  rupture  dans  LS-DYNA.  II  est  recommande  de  developper  une  methode 
‘ad  hoc’  pour  predire  l’endommagement  et  la  rupture  dans  un  modele  a  grande  echelle.  Cette 
methode  devrait  inclure  un  modele  d’endommagement,  un  critere  de  rupture,  ainsi  qu’une  etude 
de  maillage.  Des  donnees  experimentales  devraient  aussi  etre  utilisees  pour  valider  le  modele. 


ii 


DRDC  Valcartier  TM  2012-222 


Executive  summary 


Internal  Blast  in  a  Compartment-type  Vessel:  Part  1:  Finite 
Element  Modeling  Investigation 

Genevieve  Toussaint;  Claude  Fortier;  Stephane  Dumas;  DRDC  ValcartierTM 
2012-222;  Defence  R&D  Canada  -  Valcartier;  November  2012. 

Introduction  or  background:  DRDC  Atlantic  requested  DRDC  Valcartier  to  evaluate  the 
damage  on  ship  structures  such  as  bulkheads  and  decks  produced  by  a  charge  detonated  in  a 
compartment-type  vessel  by  conducting  parametric  simulations  (analytical  or  using  finite 
element)  on  charge  size,  location,  etc.  Blast  propagation  to  neighbouring  compartments  was  also 
of  interest.  After  discussions  on  the  possible  solutions  and  considering  the  short  term  deadline,  it 
was  decided  to  investigate  the  feasibility  of  using  LS-DYNA  to  perform  the  finite  element 
structural  analysis  to  study  the  structural  damage  of  a  typical  compartment-type  vessel.  The  work 
was  divided  into  two  main  activities:  first,  various  blast  loading  and  propagation  methods  usable 
in  LS-DYNA  were  investigated;  second,  a  discussion  on  the  material  damage  and  failure  models 
that  could  be  used  to  predict  the  damage  on  the  compartment  structure  is  presented. 

Results  and  significance:  This  work  has  demonstrated  the  limitations  of  actual  CFD/FE  models 
to  simulate  an  accurate  internal  blast,  including  shock  reflections  and  quasi  static  pressure,  in  a 
large  structure  and  the  need  to  address  this  gap  by  using  alternative  methods.  It  also  demonstrated 
the  need  to  develop  a  method  to  predict  the  damage  and  failure  caused  to  the  compartment 
structure. 

Future  plans:  Current  blast  propagation  methods  are  probably  adequate  for  modeling  the  quasi¬ 
static  phase  and  effects  on  neighbouring  compartments  as  long  as  the  panel  rupture  is  modeled 
realistically.  However,  on  a  more  long  term  basis,  the  development  and  validation  of  a  raytracer 
and  its  coupling  with  LS-DYNA  will  be  addressed  for  modeling  the  shock  loading  in  the  first 
compartment. 
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Internal  Blast  in  a  Compartment-type  Vessel:  Part  1:  Finite 
Element  Modeling  Investigation 

Genevieve  Toussaint;  Claude  Fortier;  Stephane  Dumas  ;  DRDC  ValcartierTM 

2012-222  ;  R  &  D  pour  la  defense  Canada  -  Valcartier;  novembre  2012. 

Introduction  ou  contexte  :  DRDC  Atlantique  a  demande  a  DRDC  Valcartier  d’evaluer  le 
dommage  cause  a  des  structures  navales  telles  que  des  cloisons  et  ponts  produites  par  une 
explosion  dans  un  compartiment  de  bateau  typique  en  effectuant  des  simulations  parametriques 
(analytiques  ou  par  elements  finis)  sur  la  dimension  de  la  charge,  sa  localisation,  etc.  La 
propagation  des  effets  de  souffle  aux  compartiments  adjacents  presentait  aussi  un  interet.  Apres 
discussion  sur  les  solutions  possibles  et  considerant  le  court  delai  pour  effectuer  le  travail,  il  a  ete 
decide  d’investiguer  la  faisabilite  d’utiliser  LS-DYNA  pour  realiser  des  analyses  par  elements 
finis  visant  a  etudier  l’endommagement  structurel  d’un  compartiment  typique  de  bateau.  Le 
travail  a  ete  divise  en  deux  activites  principals  :  en  premier,  differentes  methodes  disponibles 
dans  LS-DYNA  pour  generer  et  propager  le  souffle  d’explosion  ont  ete  etudiees;  en  deuxieme, 
une  discussion  sur  les  modeles  d’endommagement  materiel  et  de  rupture  qui  pourraient  etre 
utilises  pour  predire  l’endommagement  et  la  rupture  des  structures  est  presentee. 

Resultats  et  Importance:  Ce  travail  a  demontre  les  limitations  des  modeles  EF  actuels  pour 
simuler  adequatement  le  souffle  d’explosion,  incluant  le  choc  et  ses  reflexions,  a  l’interieur  d’une 
large  structure  et  le  besoin  de  remedier  a  cette  lacune  en  utilisant  des  methodes  alternatives.  Ce 
travail  a  aussi  demontre  le  besoin  de  developper  une  methode  pour  predire  l’endommagement  et 
la  rupture  causee  a  la  structure  du  compartiment. 

Perspectives  :  Les  methodes  actuelles  semblent  adequates  pour  modeliser  la  phase  de  pression 
quasi-statique  et  ses  effets  sur  les  compartiments  voisins,  en  autant  qu’on  dispose  de  modeles  de 
rupture  de  panneaux  realistes.  Cependant,  a  long  terme,  le  developpement  et  la  validation  d’un 
‘raytracer’  et  son  couplage  avec  LS-DYNA  seront  etudies  pour  modeliser  le  choc  dans  le 
compartiment  initial. 
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1  Introduction 


DRDC  Atlantic  requested  support  from  DRDC  Valcartier  concerning  the  Canadian  Surface 
Combatant  (CSC)  project.  The  objective  of  the  short-term1  task  to  which  this  paper  contributes 
was  to  compare  and  quantify  the  strength  and  stability  implications  of  damage  to  a  CSC-type 
concept  ship,  with  the  idea  of  highlighting  potential  design  problems.  After  discussion  between 
DRDC  Valcartier  and  DRDC  Atlantic,  internal  blast  was  identified  as  the  driving  threat  to  be 
considered. 

It  was  decided  to  conduct  parametric  simulations  (analytical  or  using  finite  element)  on  charge 
size,  location,  etc.,  to  evaluate  the  damage  on  structures  such  as  bulkheads  and  decks.  The  results 
obtained  could  be  used  to  produce  or  modify  requirements  for  the  System  Requirements 
Document  (SRD)  [Pegg  (2012)].  After  discussions  on  the  possible  solutions,  it  was  decided  to 
investigate  the  feasibility  of  using  LS-DYNA  to  perform  finite  element  structural  analyses 
studying  the  structural  damage  of  a  typical  compartment.  To  do  that,  the  work  was  divided  in  two 
main  activities:  the  first  one  studied  various  methods  used  to  generate  and  propagate  the  blast 
(including  reflections  if  necessary)  and  the  other  one  discussed  the  material  damage  and  failure 
models  in  LS-DYNA  that  could  be  used  to  predict  the  structural  behaviour  of  the  compartment 
walls,  ceiling  and  floor;  which  included  damage/erosion  on  the  elements  and  possible  venting  to 
other  compartments. 

This  work  was  performed  between  May  and  October  2012  under  task  CSC  B32  ‘Residual 
strength  and  stability  in  a  damaged  condition’. 


1  To  be  completed  in  2012 
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2  Blast  modeling 


The  first  sub-section  presents  several  methods  that  are  used  typically  in  the  DRDC  Valcartier 
Weapons  Effects  and  Protection  (WEP)  section  to  simulate  internal  and  external  blasts.  The 
second  sub-section  presents  various  finite  element  modeling  studies  that  were  performed  to 
simulate  an  internal  blast. 


2.1  Loading  methods 

Without  even  talking  about  blast  propagation  between  compartments,  simulating  the  effect  of  an 
internal  blast  in  a  single  compartment  using  the  finite  element  method  is  challenging  because  it 
implies: 

1)  modeling  multiple  shock  reflections  in  the  initial  compartment  with  sufficient  resolution 
to  represent  realistically  a  very  thin  shock  surface. 

2)  modeling  an  extensive  number  of  structural  elements  over  a  compartment  having 
dimensions  reaching  several  meters. 

Actually,  DRDC  Valcartier  uses  several  methods  to  model  land  mine  and  air  blast:  CON  WEP 
model  [Elyde  (1988),  Randers-Pehrson  and  Bannister  (1997)];  Westine  model  [Westine  (1985)]; 
Pressure -based  mine  loading  model  [Toussaint  and  Bouamoul  (2010)];  Chinook  code  [Martec 
(2005)];  arbitrary  lagrangian  eulerian  (ALE)  and  smoothed  particle  hydrodynamic  methods  (SPH) 
[Toussaint  and  Bouamoul  (2010),  LS-DYNA  Theory  Manual  (2006)].  Since  Westine  and  the  in- 
house  pressure  models  are  typically  used  to  simulate  land  mine  blast,  they  are  not  suitable  for 
internal  blast.  Also,  according  to  [Toussaint  and  Bouamoul  (2010)],  the  ALE  method  predicted 
better  the  expansion  of  air  blast  than  SPH,  so  SPH  was  not  considered.  Therefore,  to  simulate 
internal  blast  two  methods  were  considered  in  this  work:  CONWEP  and  ALE.  However,  two 
other  approaches  were  shortly  described:  Chinook  and  raytracer. 

2.1.1  CONWEP 

According  to  Schwer  (2010-A),  Kingery  and  Bulmash  (1984)  ‘ parameterized  an  extensive 
collection  of  air  blast  experimental  data  using  two  fundamental  air  blast  principals:  TNT 
Equivalence  of  difference  explosives  and  Cube  Root  Scaling  of  range,  impulse,  and  time  with 
explosive  weight’.  The  Friedlander  equation  combined  with  the  parameterization  of  these 
experimental  data,  are  the  founding  principles  of  CONWEP.  The  air  blast  portion  of  CONWEP 
was  implemented  in  LS-DYNA  by  Rahnders-Pehrson  and  Bannister  (1997)  as  an  air  blast 
function  called  *load_blast  that  treats  free-air  burst,  surface  burst  and  the  more  recent  version 
*load_blast_enhanced  (LBE),  can  also  treats  moving  free  air  burst  and  height  of  burst  type 
problem  [LS-DYNA  Keyword’s  User’s  Manual  (2012)].  One  limitation  of  the  LBE  technique  is 
that  it  can  apply  a  pressure  load  on  segments  but  the  incident  and  reflected  waves  cannot  interact. 
For  example,  LBE  cannot  account  for  reflections  produce  in  the  interior  comer  of  a  room. 
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In  CONWEP,  the  range  of  applicability  for  the  positive  phase  duration  of  a  spherical  free  air 
explosion  is  2.7  to  750  charge  radii.  That’s  why  in  simulation  runs  the  minimum  range  of 
applicability  for  CONWEP  is  usually  considered  to  be  three  times  the  charge  radii.  However, 
Schwer  (2010-B)  showed  that  this  number  was  not  strictly  correct  (see  Annex  A).  Schwer  (2010) 
presented  a  comparison  between  the  pressure  and  impulse  of  a  TNT  free  air  burst  of  4.5  kg  and 
454.5  kg  given  by  Swizdak  (1975)  and  that  predicted  by  CONWEP.  The  comparisons  between 
them  showed  that  the  maximum  incident  pressures  agreed  quite  well;  however,  CONWEP 
overpredicted  the  impulses  at  close  range.  Therefore,  a  new  range  of  applicability  of  7.57  should 
be  used  in  this  work  to  compare  the  impulse  data  with  CONWEP’s  predictions  (see  Annex  A)  for 
a  100  kg  spherical  charge  of  TNT. 

Therefore,  in  this  work,  the  *load_blast  card  in  LS-DYNA  should  not  be  used  to  load  elements  at 
a  distance  closer  than  1.852  m  for  a  100  kg  TNT  charge.  As  well,  when  we  compare  numerical 
simulation  predictions  with  CONWEP’s  impulse,  this  should  be  done  only  for  distances  further 
than  1.852  m  but  we  can  compare  the  maximum  incident  pressure  from  CONWEP  from  0.734  m 
(charge  radii  of  3). 

The  main  advantage  of  CONWEP  is  that  it  is  both  fast  (no  CFD  required)  and  reliable  when  used 
within  the  calibration  range.  The  main  disadvantage  is  that  it  is  limited  to  simple  interactions  (no 
shock  reflections). 

2.1.2  ALE 

The  ALE  formulation  allows  modeling  the  blast  as  a  fluid.  One  advantage  of  this  formulation  is 
that  it  can  model  the  interaction  of  incident  and  reflected  waves.  However,  what  limits  the  use  of 
this  model  for  large  scale  is  that  each  element  must  be  very  small  to  capture  the  shock  pressure. 
For  example,  Dong  et  al.  (2009)  modeled  air  with  element  mesh  size  smaller  than  1  mm.  Slavik 
(2009)  showed  that  an  ALE  mesh  with  element  size  fixed  to  4  mm  seemed  appropriate.  Toussaint 
and  Bouamoul  (2010)  showed  that  air  element  size  fixed  to  40  mm  seemed  adequate  to  reproduce 
the  experimental  velocity  profile  of  a  structure  resulting  from  an  air  blast.  In  the  actual  work,  the 
element  size  for  the  ALE  domain  is  one  of  the  most  difficult  parameters  to  overcome.  This  is  due 
to  the  size  of  the  compartment  (the  ALE  air  model  is  defined  by  a  box  of  dimensions:  8.1m  x 
4.6m  x  1 1.8  m).  Because  the  compartment  is  not  symmetric,  a  quarter  or  half  model  could  not  be 
used. 

2.1.3  Chinook 

Chinook23  is  a  CFD  code  designed  specifically  for  modeling  fast  blast  loadings  (shocks), 
accounting  for  phenomena  such  as  focusing,  channelling,  clearing  and  sheltering  effects.  It 
models  blast  pressure  loads  on  different  platforms  such  as  vehicles  and  personnel.  It  also  models 
multiple  materials  (explosives,  solids,  water,  gases,  and  mixtures)  and  it  models  blast  and  ejecta 
resulting  from  land  mine  explosion.  This  code  seems  well  suited  for  internal  blast  modeling, 
when  close  threat  modeling  is  required  or  for  a  complex  scenario  modeling  and  was  validated 
with  experimental  data.  Martec  has  previously  used  Chinook  for  CPF  vulnerability  to  IEDs 


2  Chinook  Software  was  developed  by  Martec  Limited  in  partnership  with  DRDC. 

3  www.martec.com 
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[Martec  (2005)].  Chinook  is  seen  as  more  efficient  than  LS-DYNA’s  internal  CFD-like  module 
(ALE,  SPH,  CONWEP)  for  modeling  shock  effects  and  may  be  coupled  with  LS-DYNA  for 
modeling  fluid-structure  interactions.  The  2D  code  seems  to  work  with  reasonable  run  times; 
however,  the  execution  time  for  the  3D  code  is  very  long.  Therefore,  we  may  encounter  the  same 
limitations  as  in  the  3D  ALE  model  concerning  the  small  cell  size  required  to  model  accurate 
maximum  incident  pressure. 

Since  coupling  with  Chinook  has  not  been  fully  experimented  yet  by  DRDC  Valcartier  and  due  to 
the  timeline  given  to  realise  the  work,  Chinook  was  only  explored.  In  order  to  get  started  using 
Chinook  and  for  future  needs,  a  training  session  was  organised  in  September  2012.  This  included 
exploration  of  the  fluid-structure  coupling. 

2.1.4  Raytracer 

The  raytracer  is  a  method  using  virtual  sources  to  assess  the  pressure  at  a  location  (i.e.  a  gage). 
The  version  presently  experimented  assumes  a  rectangular  room  with  perpendicular  walls.  The 
algorithm  can  support  any  number  of  explosives  and  gages.  The  resultant  pressures  recorded  is 
the  combined  pressure  histories  from  the  incident  shockwaves  and  reflections  from  the  walls  and 
from  all  explosives  (several  at  a  time  may  be  modeled).  Another  advantage  is  that  it  avoids  the 
modeling  of  air  cells  so  that  loading  the  walls  is  much  faster  than  using  ALE  or  Chinook.  This 
method  is  still  under  development  and  coupling  with  LS-DYNA  has  to  been  done  yet.  However, 
it  is  a  promising  tool  and  should  allow  efficient  and  reasonably  accurate  modeling  of  an  internal 
blast  in  a  compartment-type  vessel. 

2.2  Finite  Element  Modeling 

This  sub-section  presents  four  studies  that  were  realised  to  investigate  the  feasibility  of  using  LS- 
DYNA  to  simulate  accurately  the  propagation  of  an  internal  blast  considering  the  size  of  a  typical 
compartment-type  vessel.  The  first  study  presents  the  simulation  of  a  blast  propagating  in  a 
simple  multi-material  ALE  FE  model  of  a  rectangular  box  meshed  with  elements  of  200  mm.  In 
the  second  study,  this  FE  model  was  refined  and  three  new  FE  models  were  generated  with 
element  sizes  of  100  mm,  75  mm  and  50  mm.  The  third  study  presents  the  combination  of  the 
*load_b!ast_enhanced  (LBE)  keyword  and  the  multi-material  ALE  approach  to  generate  the 
blast.  The  fourth  study  presents  the  propagation  of  the  blast  in  a  spherical  domain.  Finally, 
alternative  approaches  to  generate  and  propagate  the  blast  are  presented. 

Note  that  in  principle,  cells  should  be  taken  much  smaller  than  those  above  to  capture  shock 
features.  However,  higher  resolutions  would  be  almost  unmanageable,  especially  if  extended  to 
several  compartments.  It  is  then  of  interest  to  check  if  convincing  loadings  could  be  modeled  with 
the  above  resolutions.  Given  the  relatively  slow  panel  response  times  compared  to  shock 
durations,  impulse  may  be  more  important  than  (very  short  duration)  peak  pressures. 

The  compartment-type  vessel,  as  shown  in  Figure  1,  is  made  of  shell  elements  and  beam  sections. 
Notice  that  the  compartment  is  not  symmetric.  Notice  also  that,  the  final  objective  is  to  be  able  to 
use  the  finite  element  model  to  conduct  parametric  simulations  on  charge  size,  location  (ex. :  top 
corner,  bottom  corner,  mid-wall...),  etc.,  to  evaluate  the  damage  on  structures  such  as  bulkheads 
and  decks.  Therefore,  the  approach  chosen  should  allow  the  positioning  of  the  charges  easily 
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everywhere  in  the  compartment  without  the  need  to  regenerate  the  FE  model  for  each  possible 
scenario. 


Figure  1:  Compartment-type  vessel 


2.2.1  First  study 

The  first  FE  model  created  in  FEMAP4 5  was  a  rectangular  box  with  dimensions  of  5.8  m  x  9.6  m  x 
13.2  m.  A  total  of  91,872  elements  were  generated  with  element  size  of  200  mm.  A  sphere 
representing  a  mass  of  100  kg  TNT  was  defined  using  the  *lnitial_volume __ fraction  geometry 
card  and  positioned  in  the  middle  of  the  box.  The  mesh  was  uniform  because  the  FE  model  had  to 
be  valid  for  any  charge  locations  in  the  box. 


Figure  2:  ALE  rectangular  mesh5 with  element  size  of 200  mm 


4  FEMAP  was  used  to  generate  all  the  FE  models  of  this  work. 

5  The  units  for  all  the  FE  models  and  results  in  this  work  are:  g,  mm,  ms,  MPa,  N. 
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A  pressure  of  101  kPa  was  applied  on  the  outside  boundaries  of  the  box  and  inside  the 
compartment.  A  ^boundary  non  reflecting  (BNR)  card  was  specified  to  avoid  reflection  waves 
to  re-enter  in  the  model.  The  air  was  modeled  with  the  *mat_null  material  model  using  a  linear 
polynomial  equation  of  state.  Properties  are  provided  in  Table  1  and  Table  2.  The  explosive  was 
modeled  using  *  mat _h ighexp l osi veb urn  material  combined  with  the  Jones  Wilkins  Lee  (JWL) 
equation  of  state.  Properties  are  provided  in  Table  3  and  Table  4. 

Table  1:  *Mat  null  material  parameters  [Toussaint  and  Bouamoul  (2010)] 


Air 

Density  (p),  g/mm3  1.293xl0'6 

Pressure  cutoff  (Pc),  MPa  0 


Table  2:  *EOS  linear  polynomial  [Toussaint  and  Bouamoul  (2010)] 

Air 

Polynomial  equation  coefficient  CO,  Cl,  C2,  C3  and  C6  0.0 
Polynomial  equation  coefficient  C4  and  C5  0.4 

Initial  internal  energy  (E0),  MPa  0.25 

Initial  relative  volume  (V0)  1.0 


Table  3:  Explosive  properties  [Dobratz  and  Crawford  (1985)] 


TNT 

Density  (p),  g/mm3 

1.630xl0'3 

Detonation  velocity  (D),  mm/ms 

6930 

Chapman-Jouget  pressure  (Pcj),  MPa 

21xl03 

Table  4:  *EOS  JWL  [Dobratz  and  Crawford  (1985)] 

JWL 

A,  MPa 

3.712xl05 

B,  MPa 

3.231 

xlO3 

R1 

4.15 

R2 

0.95 

Omega  (oo) 

0.30 

Internal  energy  density  (E0),  MPa  - 

7.0xl03 

mm3/mm3 

Initial  relative  volume  (Vo) 

1.0 
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The  finite  element  structural  code  LS-DYNA  was  used  to  simulate  the  blast  wave  propagation  in 
the  box. 

For  the  first  run,  the  center  of  detonation  was  located  in  the  middle  of  one  element  (y  axis)  as 
shown  in  Figure  3.  Also,  one  effect  of  having  a  coarse  mesh  is  that  the  explosive  is  represented  by 
a  tetrahedron  instead  of  a  sphere  made  of  many  cells. 


Figure  3:  a)  Location  of  the  center  of  detonation  b)  Explosive  representation  in  ALE  model 

When  analysing  the  results  of  the  simulation,  we  observed  that  the  results  obtained  on  the  axis  (y) 
were  different  than  the  ones  obtained  on  the  x  and  z  axis  (at  the  same  range).  To  reduce  this 
effect,  the  center  of  the  explosion  and  the  detonation  point  were  moved  from  the  initial  location 
on  the  y  axis  to  the  closest  node  location,  as  shown  in  red,  in  Figure  4.  Even  by  doing  so,  at  0°  on 
the  xy,  yz  and  xz  plane,  the  pressures  predicted  were  slightly  different  at  the  same  range  from  the 
center  of  detonation.  This  could  be  explained  by  the  fact  that  even  if  a  BNR  card  was  used,  there 
were  still  reflections  at  the  boundaries.  This  could  be  reduced  by  pushing  the  boundaries  far 
enough  from  the  center  of  detonation  but  in  this  case,  it  was  not  possible  due  to  the  size  of  the 
problem. 


Figure  4:  New’  location  of  the  center  of  detonation 


Giving  the  range  of  applicability  from  Swizdak  (1975),  for  a  100  kg  charge,  the  distance  in 
charge  radii  should  be  set  7.57  (see  Annex  A)  which  corresponds  to  1.852  m  for  the  comparison 
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between  the  impulse  predictions  from  CONWEP  and  numerical  results.  The  maximum  incident 
pressure  and  impulse  obtained  at  1.5  m,  2  m  and  3  m  are  given  in  Table  5  and  Table  6. 


Table  5:  Maximum  incident  pressure  obtained  for  a  rectangular  mesh  with  element  size  of  200 

mm  compared  to  the  CONWEP  predictions 


Distance  (m) 

Mesh  200  mm 

CONWEP 

Pressure  (MPa) 

Pressure  (MPa) 

1.5 

14.148 

8.059 

2.0 

5.07413 

5.064 

3.0 

4.626053 

2.369 

Table  6:  Impulse  obtained for  a  rectangular  mesh  with  element  size  of 200  mm  compared  to  the 

CONWEP  predictions 


Distance  (m) 

Mesh  200mm 

CONWEP 

Impulse 

(MPa«ms) 

Impulse 

(MPa*ms) 

1.5 

1.27703 

2.0 

0.836873 

0.6352 

3.0 

1.2140875 

0.7585 

Note  that  when  the  location  of  the  pressure  gauge  fell  between  two  elements,  the  mean  of  the 
pressure  of  both  elements  was  taken. 

This  table  highlights  some  limitations  of  this  model;  the  FE  model  overpredicts  the  maximum 
incident  pressure  and  impulse  compared  to  CONWEP.  In  order  to  see  if  refining  the  mesh  will 
increase  the  maximum  incident  pressure  and  total  impulse,  a  convergence  study  is  required  and 
will  be  presented  in  the  next  section. 

Lessons  were  learned  from  this  study: 

The  center  of  detonation  must  be  coincident  with  a  node  (not  in  the  middle  of  an 
element). 

Even  if  a  non-reflective  condition  is  given  on  the  outside  boundaries,  the  boundaries 
should  be  as  far  as  possible. 

Elements  must  be  smaller  than  200  mm. 

The  explosive  should  be  defined  with  many  elements  (depending  on  the  charge  radius). 
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2.2.2  Second  study 


In  the  previous  sub-section,  a  3D  rectangular  mesh  was  generated  with  a  mesh  size  of  200  mm.  It 
was  demonstrated  that  elements  with  a  mesh  size  of  200  mm  did  not  capture  accurately  the 
maximum  incident  pressure  and  impulse  at  different  locations  from  the  center  of  detonation. 
Also,  it  was  showed  that  the  size  of  the  mesh  cells  must  be  sufficiently  small  to  allow  the 
propagation  of  a  nice  spherical  blast. 

In  this  study,  a  finite  element  model  of  dimensions  8.1m  x  4.6m  x  11.8  m  was  modeled  with 
hexahedron  ALE  elements  with  mesh  sizes  of  100  mm  (439,668  elements),  75  mm  (1,034,316 
elements)  and  50  mm  (3,  517,344  elements).  Elements  were  not  smaller  than  50  mm  because  it 
generated  too  many  elements  that  couldn’t  be  handled  by  FEMAP  on  my  workstation6.  A 
boundary  non-reflecting  card  was  specified  to  limit  the  reflection  at  the  boundaries.  The  meshes 
are  shown  in  the  following  figures. 


2 

b 

4 

L 


Figure  6:  Rectangular  mesh-  75  mm 


6  We  tried  modeling  the  rectangular  mesh  with  elements  size  of  40  mm  without  success. 
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Figure  7:  Rectangular  mesh-50  mm 

The  finite  element  structural  code  LS-DYNA  v5.1.1  using  multiple  processors  was  used  to 
simulate  the  blast  wave  propagation  of  a  100  kg  TNT  charge  located  in  the  middle  of  the  box. 
Several  problems  were  encountered.  For  example,  the  timestep  of  the  d3plot  recorded  needed  to 
be  decreased  in  order  to  capture  the  peak  pressure  that  fell  between  two  measurements  in  the 
simulation.  The  center  of  detonation  had  to  be  moved  to  fit  with  a  node  location  (see  Figure  5  and 
Figure  6)  such  as  in  the  previous  study.  Also,  there  were  some  reflections  in  the  model  even  if  a 
BNR  card  was  used  and  boundaries  could  not  be  moved  further.  Figure  8  shows  an  example  of  an 
air  blast  propagation. 
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Figure  8:  Example  of  blast  propagation  in  the  FE  model  (MPa) 


A  comparison  of  the  maximum  incident  pressure  and  total  impulse  taken  at  1.5  m,  2  m  and  3  m 
obtained  for  the  different  meshes  and  the  CONWEP  predictions  are  given  in  Table  7  and  Table  8. 
Annex  B  presents  an  example  of  the  curves  obtained. 
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Table  7:  Maximum  incident  pressure  obtained  for  a  rectangular  mesh  with  element  size  of  100 
mm,  75  mm  and  50  mm  compared  to  CONWEP  results 


Distance  (m) 

Peak  Pressure  (MPa) 

CONWEP 

(MPa) 

Mesh  100  mm 

Mesh  75  mm 

Mesh  50  mm 

1.5 

1.878 

2.972 

4.357 

8.059 

2.0 

1.861 

2.325 

2.957 

5.064 

3.0 

1.195 

1.446 

1.712 

2.369 

Table  8:  Impulse  obtained for  a  rectangular  mesh  with  element  size  of  100  mm,  75  mm  and  50 

mm  compared  to  CONWEP  results 


Distance  (m) 

Impulse  (MPa«ms) 

CONWEP 

(MPa*ms) 

Mesh  100  mm 

Mesh  75  mm 

Mesh  50  mm 

1.5 

0.2997 

0.3573 

0.4260 

0.6723 

2.0 

0.3886 

0.4372 

0.4926 

0.6352 

3.0 

0.6436 

0.6687 

0.6752 

0.7585 

Note  that  when  the  location  of  the  pressure  fell  between  two  elements,  the  mean  of  the  pressure  of 
both  elements  was  taken. 

Table  7  shows  that  peak  pressure  increases  when  refining  the  mesh  and  decreases  when 
increasing  the  distance.  It  can  also  be  seen  that  a  mesh  of  50  mm  is  not  sufficient  to  capture 
adequately  the  maximum  incident  pressure  predicted  by  CONWEP.  According  to  Mahmadi  et  al. 
(2002)  and  Alia  and  Souli  (2006),  it  is  possible  to  match  accurately  the  peak  pressure  when  the 
mesh  is  fine  enough.  In  this  case,  it  would  probably  require  running  the  problem  with  element  in 
the  order  of  1  to  1 0  mm  which  seems  almost  impossible  for  a  ship  size  problem. 

Table  8  shows  that  total  impulse  increases  when  refining  the  mesh  and  increases  also  when  the 
distance  increases.  Results  obtained  are  still  not  sufficiently  close  to  the  CONWEP  predictions. 
Moreover,  CONWEP  predicts  a  reduction  of  the  total  impulse  at  2  m  that  is  not  captured  by  any 
of  the  three  models. 

To  reduce  the  multi-material  ALE  domain  (i.e.  reduce  the  number  of  elements),  another  approach 
has  to  be  chosen.  This  method  is  presented  in  the  next  section  and  consists  in  coupling  the  multi¬ 
material  ALE  model  with  a  *load_blast_enhanced  keyword  in  LS-DYNA  as  presented  in  Schwer 
(2010-B). 

Some  lessons  were  learned  from  this  study: 

Elements  smaller  than  50  mm  are  required  to  capture  peak  pressure  and  total  impulse 
produced  by  a  100  kg  TNT  detonation  at  close  range  (<  to  3  m). 
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There  is  a  need  to  increase  the  frequency  of  d3plot  data  at  a  minimum  of  0.01ms  to 
capture  the  highest  pressure. 

The  pressure  in  the  elements  along  each  principal  axes  and  at  45°  should  always  be 
compared. 

2.2.3  Third  study 

The  previous  sub-section  showed  that  the  element  mesh  must  be  smaller  than  50  mm  to  capture 
the  maximum  incident  pressure  and  impulse  data.  The  objective  of  this  third  approach  is  to  find  a 
way  to  reduce  the  number  of  ALE  elements  in  the  model  and  consequently,  to  increase  the  mesh 
size  to  get  maximum  incident  pressure  and  total  impulse  closer  to  the  CONWEP  predictions.  One 
solution  to  limit  the  number  of  elements  in  the  model,  and  at  the  same  time  minimize  the 
computational  cost,  is  to  combine  the  LBE  technique  with  an  ALE  model.  This  technique  is  well 
explained  in  Schwer  (2010-B).  The  load  produced  by  LBE  is  applied  to  the  elements  at  the  outer 
surface  of  the  eulerian  mesh.  The  interior  wave  reflections  are  then  treated  by  the  multi-material 
ALE  solver.  With  this  method,  there  is  no  need  to  model  all  the  air  elements  from  the  charge  to 
the  target,  only  a  smaller  air  box  is  required  around  each  wall/ceiling. 

In  the  multi-material  ALE  approach  coupled  to  the  *load_blast_enhanced  (LBE-ALE)  keyword 
in  LS-DYNA,  the  explosive  is  not  modeled  physically;  instead,  a  load  pressure  is  applied  to  the 
air  segments  forming  the  boundary  facing  the  location  of  the  charge.  Figure  9  shows  an  example. 
Instead  of  modeling  the  whole  domain  with  multi-material  ALE  (like  in  section  2.2),  only  a 
portion  of  the  model  is  meshed  (pink  zone)  and  the  pressure  is  applied  on  the  first  row  of 
elements  forming  the  frontier.  One  advantage  of  using  this  technique  is  to  reduce  the  domain 
(number  of  elements)  to  be  modeled  and  another  advantage  is  that  reflections  can  still  occur  in  the 
multi-material  ALE  domain.  However,  outside  the  multi-material  ALE  domain,  reflections  are 
not  accounted  for.  Also,  in  this  work,  LBE  should  not  be  applied  on  a  frontier  closer  to  3  charge 
radii  (for  maximum  incident  pressure)  and  preferably  8  charge  radii  (for  impulse)  which  means 
that  near  field  blast  and  multiple  internal  blasts  cannot  be  modeled  with  this  technique. 


Figure  9:  LBE-ALE  example 

In  this  study,  since  only  a  portion  of  the  3D  rectangular  box  presented  previously  was  used,  the 
mesh  could  be  refined.  Figure  10  shows  the  LBE  frontier  (in  purple)  for  three  meshes  size  of  100 
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mm,  50  mm  and  10  mm.  For  each  mesh,  the  ALE  model  (including  the  first  layer  of  elements 
where  the  LBE  is  applied)  was  moved  at  the  different  locations  (1.5  m,  2  m  and  3  m  from  the 
center  of  detonation)  to  get  maximum  incident  pressure  and  impulse  at  each  of  these  locations. 


Figure  10:  ALE  mesh  where  the  load  is  applied  (100  mm  -  50  mm  -  10  mm) 


The  finite  element  structural  code  LS-DYNA  v5.1.1  was  run  on  a  single  processor  to  simulate  the 
blast  wave  propagation  of  a  100  kg  TNT  charge  located  in  the  middle  of  the  box.  Note  that  the 
loading  generated  on  these  segments  (layer)  did  not  depend  on  the  number  of  elements  in  the 
back.  Table  9  and  Table  10  show  the  results  of  the  peak  pressure  and  impulse  measured  at  1.5  m, 
2  m  and  3  m  for  each  LBE-ALE  mesh  (measurements  taken  perpendicularly  to  the  charge).  The 
results  obtained  are  compared  to  the  CONWEP  predictions. 


Table  9:  Maximum  incident  pressure  obtained  for  the  LBE-ALE  model  with  mesh  sizes  of  100 
mm,  50  mm  and  10  mm  compared  to  the  CONWEP  predictions 


Distance  (m) 

Peak  Pressure  (MPa) 

CONWEP 

(MPa) 

Mesh  100  mm 

Mesh  50  mm 

Mesh  10  mm 

1.5 

7.648 

7.617 

6.823 

8.059 

2.0 

4.535 

4.491 

4.322 

5.064 

3.0 

2.328 

2.326 

2.175 

2.369 

Table  10:  Impulse  obtained for  the  LBE-ALE  model  with  mesh  sizes  of  100  mm,  50  mm  and  10 

mm  compared  to  the  CONWEP  predictions 


Distance  (m) 

Impulse  (MPa*ms) 

CONWEP 

(MPa»ms) 

Mesh  100  mm 

Mesh  50  mm 

Mesh  10  mm 

1.5 

0.6857 

0.6714 

0.6289 

0.6723 

2.0 

0.6095 

0.6027 

0.5772 

0.6352 

3.0 

0.7525 

0.7440 

0.7293 

0.7585 
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Table  9  shows  that  the  peak  pressure  decreases  when  refining  the  mesh.  Table  10  shows  that  the 
impulse  decreases  when  refining  the  mesh.  We  expected  the  opposite  behaviour. 

For  the  same  element  size,  the  FE  model  results  are  much  closer  to  the  CONWEP  predictions 
when  using  the  LBE-ALE  model  than  using  only  a  multi-material  ALE  model. 

Both  tables  showed  that  the  best  correlation  between  the  FE  model  and  the  CONWEP  predictions 
is  for  a  mesh  with  element  size  of  100  mm.  Therefore,  with  element  size  of  100  mm,  the  LBE- 
ALE  captured  adequately  the  initial  loading.  This  approach  could  be  used  to  model  the  effect  of  a 
direct  loading  only  (no  reflections)  on  the  compartment  walls.  In  this  work,  it  will  not  be  useful  as 
we  want  to  model  the  effect  of  multiple  reflections  on  the  compartment  structure. 

2.2.4  Fourth  study 

In  the  previous  studies,  a  rectangular  mesh  domain  was  used  to  propagate  the  blast.  The  main 
disadvantage  of  a  rectangular  grid  is  that  even  at  the  same  distance  from  the  center  of  detonation, 
the  pressure  predicted  varied  in  all  directions.  Schwer  et  al.  (2008)  compared  the  maximum 
impulse  in  the  xy  plane  at  0°,  22.5°,  45°,  67.5°  and  90°  and  showed  large  variations  between  the 
results  when  using  a  rectangular  mesh  compared  to  a  spherical  one.  It  is  possible  to  reduce  this 
directional  dependence  by  refining  the  mesh  or  to  eliminate  the  problem  by  using  a  spherical 
coordinate  system. 

In  this  section  we  generated  a  spherical  domain.  One  advantage  of  using  a  spherical  domain  is 
that  it  requires  a  lot  less  element  than  using  a  rectangular  domain  so  the  elements  near  the 
detonation  can  be  smaller.  Another  advantage  is  that  there  is  no  directional  dependence  of  the 
mesh.  However,  once  the  blast  has  reached  the  structure  the  wave  reflected  will  not  be  in  the 
direction  of  the  spherical  mesh  and  therefore  the  reflections  might  not  be  accurate7.  Also,  using 
this  mesh,  only  a  centered  detonation  in  the  compartment  is  possible,  because  an  off  centered 
detonation  (in  the  corner  for  example)  would  require  moving  the  mesh  in  the  corner  and  would 
require  generating  more  elements  to  the  domain  to  enclose  all  the  compartment  structure.  This 
would  add  too  many  elements  to  the  mesh. 

The  finite  element  model  of  the  explosive  modeled  was  generated  using  a  small  cube  box  in  the 
center  of  a  sphere  with  a  radius  of  245  mm.  For  generating  the  FE  model  of  the  air,  the  elements 
from  the  surface  of  the  explosive  were  extruded  using  300  elements  on  a  4.755  m  radius  with  a 
thickness  bias  of  5.  Three  element  sizes  of  10  mm,  5  mm  and  4  mm  for  the  initial  cubic  box  were 
generated.  The  initial  mesh  of  4  mm  had  a  higher  growth  factor  than  the  initial  mesh  of  5mm  in 
order  to  keep  the  number  of  elements  acceptable  (5  mm:  3,222,400  elements  and  4  mm: 
3,290,000  elements). 

The  central  5  mm  mesh  is  shown  in  Figure  11.  Due  to  the  number  of  elements  generated,  the 
radius  of  the  domain  had  to  be  set  to  a  maximum  of  5  m  which  allowed  to  include  only  the 
nearest  walls  of  the  compartment,  not  the  whole  compartment.  A  boundary  non-reflecting  card 
was  specified  to  limit  the  reflection  at  the  boundaries. 


7  It  might  be  possible  to  get  around  the  problem  by  using  an  adaptative  mesh,  but  it  was  not  tested  in  this 
work. 
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Figure  11:  Initial  mesh:  5  mm 


The  finite  element  structural  code  LS-DYNA  v5.1.1  was  used  to  simulate  the  blast  wave 
propagation  of  a  100  kg  TNT  charge  located  in  the  middle  of  the  sphere.  Many  FE  analyses  were 
performed  to  study  the  effect  of  parameters  such  as:  element  size,  growth  factor,  equation  of  state 
for  the  explosive  and  to  study  the  influence  of  using  method  2  versus  method  3  in  the 
*control_ale  card.  All  the  results  obtained  from  these  FE  analyses  are  not  presented  in  this  work, 
only  the  final  results  for  the  mesh  size  of  5  mm  using  the  JWL  equation  of  state  and  method  2  in 
the  *control_ale  card  are  presented.  These  results  are  given  in  Table  1 1  and  Table  12. 

Table  11:  Maximum  incident  pressure  obtained  for  a  spherical  model  with  mesh  sizes  of  5  mm 

compared  to  CONWEP  predictions 


Distance  (m) 

Mesh  5  mm 

CONWEP 

Pressure 

(MPa) 

Pressure 

(MPa) 

1.5 

5.9797 

8.059 

2.0 

4.1897 

5.064 

3.0 

2.0650 

2.369 

Table  12:  Impulse  obtained for  a  spherical  model  with  mesh  sizes  of  5  mm  compared  to 

CONWEP  predictions 


Distance  (m) 

Mesh  5mm 

CONWEP 

Impulse 

(MPa*ms) 

Impulse 

(MPa*ms) 

1.5 

0.6829 

0.6723 

2.0 

0.6864 

0.6352 

3.0 

0.6246 

0.7585 

Table  1 1  shows  that  maximum  incident  pressures  predicted  by  the  spherical  multi-material  ALE 
mesh  are  26%  lower  at  1.5  m,  17%  lower  at  2  m  and  13%  at  3  m  than  the  CONWEP  predictions. 
Flowever,  when  these  predictions  are  compared  to  the  predictions  of  the  multi-material  ALE 
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rectangular  mesh  with  the  smallest  mesh  (Table  7),  the  spherical  mesh  provides  a  better 
agreement  with  the  CONWEP  predictions. 

Table  12  shows  that  impulse  predicted  by  the  spherical  multi-material  ALE  mesh  (after  1.852  m) 
are  8%  higher  at  2  m  and  18%  lower  at  3  m  than  the  CONWEP  predictions. 

One  advantage  of  the  spherical  mesh  is  that  the  maximum  incident  pressure  and  impulse 
predicted  by  the  spherical  domain  compared  to  the  CONWEP  predictions  were  closer  than  the 
ones  obtained  with  a  rectangular  domain  having  approximately  the  same  number  of  elements. 
Another  advantage  of  the  spherical  mesh  over  the  rectangular  mesh  was  that  there  is  no 
directional  dependence  of  the  mesh.  However;  in  a  spherical  mesh,  the  loading  resulting  from 
multiple  reflections  would  not  be  accurate. 

2.2.5  Alternative  approaches 

The  approaches  taken  to  simulate  the  blast  wave  propagation  in  a  compartment  type  vessel  using 
LS-DYNA  were  not  very  successful.  Other  methods  could  have  been  tested.  However,  due  to  the 
short  time  given  to  realise  this  work,  it  was  not  possible  to  investigate  all  of  them.  Here  are  a  few 
suggestions  that  could  be  tested  in  the  future: 

1.  Divide  the  rectangular  mesh  into  10  parts,  each  part  having  its  own  input  file.  The 
element  size  in  each  part  could  be  around  25  mm.  In  the  main  file,  use  the  * include  card 
to  include  all  the  parts  in  one  model  and  use  the  *nodes_merge_set  option  in  LS-DYNA 
to  join  all  the  nodes  between  the  parts.  Instead  of  being  limited  by  FEMAP,  it  is  the 
execution  time  that  will  be  limited  by  the  number  of  elements. 

2.  One  possible  solution  to  reduce  the  number  of  elements  in  the  model  is  demonstrated  in 
Aquelet  and  Souli  (2008).  The  technique  is  to  map  axi-symmetric  blast  wave  propagation 
from  2D  ALE  mesh  to  3D  ALE  mesh  and  then,  to  run  the  ALE  3D  problem  with  LS- 
DYNA  code.  Zakrisson  et  al.  (2011)  used  that  technique  to  simulate  a  near-field  0.75  kg 
charge  weight  detonated  in  free  air.  A  convergence  study  was  performed  to  validate  the 
choice  of  the  2D  and  3D  mesh.  They  ended  using  element  side  length  of  0.5  mm  for  the 
2D  mesh  mapped  to  4  mm  element  length  (using  a  biaised  distribution)  for  the  3D  mesh. 
Mapping  from  2D  to  3D  ALE  mesh  appears  feasible  in  the  actual  work;  it  would  allow 
modeling  a  very  fine  mesh  in  two  dimensions  before  the  mapping  to  3D.  Also,  it  would 
require  two  convergence  studies  instead  of  one.  Since  the  ceiling  and  floor  are  located 
approximately  2  m  from  the  center,  the  model  would  be  useful  for  a  centered  detonation 
between  ceiling  and  floor  only  and  not  too  close  from  the  walls.  Of  course,  the  remaining 
parts  of  the  domain  would  have  to  be  modeled  with  3D  elements.  This  would  also  solve 
only  part  of  the  problem  as  subsequent  reflections  would  not  be  accurate. 

3.  Map  2D  mesh  to  a  3D  mesh  using  Chinook.  The  same  limitations  as  stated  in  the 
previous  point  form  would  apply.  A  full  3D  Chinook  solution,  would  very  likely  require 
resolutions  too  high  to  be  manageable  on  such  a  large  scale  (i.e.  same  problems  as  with 
ALE).  This  was  discussed  with  Martec  specialists  that  concluded  that  modeling  multiple 
reflections  until  the  onset  of  quasi-static  pressure  was  not  realistic. 
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4.  Finally,  the  most  promising  approach  seems  to  be  to  use  a  raytracer  to  provide  the 
loading  inside  the  compartment  in  a  very  fast  but  reasonably  accurate  manner,  accounting 
for  multiple  reflections8. 


2.3  Conclusion 

In  this  section,  we  investigated  the  feasibility  of  using  LS-DYNA  to  simulate  accurately  the 
propagation  of  an  internal  blast  in  a  large  domain  (equivalent  to  the  size  of  a  typical 
compartment-type  vessel).  Several  approaches  were  presented  and  demonstrated  that  there  are 
many  constraints  that  need  to  be  addressed  before  getting  accurate  predictions.  A  full  3D 
CFD/ALE  modeling  of  the  loading,  including  multiple  reflections,  with  a  resolution  sufficient  to 
represent  shock  effects  accurately  everywhere  in  the  compartment,  is  almost  certainly  unrealistic. 
However,  alternative  approaches  were  suggested. 


8  At  this  time,  work  was  started  to  program  the  raytracer,  but  some  validation  has  to  be  done  and  the 
coupling  with  LS-DYNA  has  not  been  implemented  yet.  This  is  part  of  future  work  that  could  be  addressed 
in  future  CSC  stemming  from  B-31[Peg  et  al.  (2012)]. 
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3  Damage  and  failure  modeling 


This  section  highlights  the  most  important  aspects  to  consider  when  modeling  damage/failure  of 
elements  in  a  FE  element  model.  First,  we  looked  at  the  way  LS-DYNA  achieves  erosion  in  its 
material  models.  Then,  we  present  only  three  studies  on  the  modeling  of  damage  of  bulkheads 
and  decks  as  there  was  very  limited  amount  of  information  available  in  the  open  literature  on  this 
subject.  These  studies  could  be  very  useful  if  this  project  is  pursued. 

3.1  Erosion/failure  in  LS-DYNA 

According  to  Schwer  (2009),  one  possible  solution  to  include  erosion  in  a  finite  element  analysis 
is  to  follow  these  three  steps: 

1 .  Select  an  ad  hoc  erosion  criterion 

2.  Select  a  critical  value  for  the  erosion  criterion 

3.  Perform  a  mesh  study:  the  results  should  be  converging  and  the  error  due  to  discretization 
should  be  small  enough. 

Since  the  selection  is  ad  hoc,  the  erosion  criterion  should  be  calibrated  with  experimental  data  or 
at  least,  a  mesh  refinement  analysis  should  be  performed  on  the  model. 

The  most  common  criterion  used  in  LS-DYNA  is  the  effective  plastic  strain.  Since  the  value  for 
erosion  in  the  element  is  the  same  in  tension  and  in  compression,  it  is  not  recommended  to  choose 
this  type  of  criterion  unless  experimental  data  are  available  to  validate  the  criterion.  In  Martec 
(2005)  report  for  example,  material  properties  were  modified  to  include  plasticity  and  failure 
strain. 

In  LS-DYNA  when  using  the  *mat  add  erosion  card,  it  is  possible  to  define  multiple 
independent  criterion  based  on,  for  example,  the  maximum  or  minimum  pressure,  the  maximum 
or  minimum  principal  strain,  the  maximum  shear  strain,  etc.  It  is  also  possible  to  include  a 
GISSMO  (Generalized  Incremental  Stress  State  dependant  damage  Model)  damage  accumulation 
model  including  softening  and  failure  or  the  user  may  invoke  an  arbitrary  number  of  damage 
initiation  and  evolution  criteria  [LS-DYNA  Keyword  User’s  Manual  Version  971  volume  I  and  II 
(2012)]. 

When  damage  or  failure  is  included  in  a  simulation,  the  mesh  topology  is  an  important  aspect 
because,  damage  and  failure  are  often  initiating  in  the  smallest  elements  of  a  mesh.  Therefore, 
doing  a  mesh  refinement  study  to  establish  the  failure  criterion  on  the  compartment  is  required.  A 
mesh  as  uniform  as  possible  should  be  used  when  predicting  erosion,  otherwise  the  smallest 
elements  should  be  modeled  in  the  area  where  damage  and  failure  seems  probable.  By  changing 
the  size  of  element  at  different  locations,  we  might  change  the  predicted  results.  Experimental 
data  seems  required  to  validate  the  failure  model. 
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3.2  Studies  on  the  modeling  of  damage  of  bulkheads  and 
decks 

Iwahashi  et  al.  (1998)  presented  a  comparative  study  of  the  local  stresses  developed  at  the 
intersection  structure  of  longitudinal  stiffeners  and  a  transverse  bulkhead  between  eight  shell 
element  models  and  two  solid  models.  The  author  recommended  using  3D  solid  elements  to  get  a 
precise  prediction  of  local  stresses  at  weld  connections.  If  shell  elements  were  to  be  used,  a 
reliable  and  practical  procedure  has  to  be  developed  to  allow  a  proper  prediction  of  the  local 
stresses. 

Raymond  (2001)  presented  in  his  master  of  engineering  the  tools  needed  for  the  formation  of 
optimised  X-80  steel  blast  tolerant  transverse  bulkheads.  This  thesis  included  all  the  procedure 
followed  to  get  to  this  optimisation.  Amongst  other  things,  the  definition  of  design  criteria  from 
the  worst  case  operational  requirements,  high  strain  rate  data  and  analysis  to  determine  the 
constants  and  constitutive  models  to  be  used  in  the  FE  model,  development  and  evaluation  of  a 
FE  modelling  technique  were  presented.  The  thesis  presented  some  factors  related  to  the  design 
constraints  (stress  waves  and  strain  rates  and  prediction  of  rupture)  and  an  investigation  on 
bulkhead  components  (joint  and  stiffener  structures). 

Tyler-Stree  and  Luyten  (2009)  have  used  the  dedicated  element  methodology  (DEM)  developed 
by  Dillingh  (2003)  specifically  to  predict  ductile  failure  of  welded  bulkhead-deck  connections 
due  to  an  internal  blast.  The  authors  reported  that  due  to  the  size  of  the  structure  involved  (ship), 
it  was  not  possible  to  use  a  fine  element  mesh  to  model  the  structure.  Therefore,  to  get  a 
reasonable  amount  of  energy  dissipation  and  a  more  realistic  elongation  before  failure,  the  failure 
criteria  should  be  adjusted  to  larger  element  size  (i.e.  in  the  order  of  typical  shell  meshes  used  for 
ship  analyses). 

3.3  Conclusion 

Due  to  time  constraint,  it  was  not  possible  to  perform  FE  analysis  including  damage  and  failure  to 
evaluate  the  damage  on  structures  (such  as  bulkheads  and  deck)  in  the  compartment-type  vessel. 
Flowever,  in  the  future,  a  method  should  be  developed  to  predict  the  damage  and  failure  in  large 
scale  model. 

This  method  could  be  to  use  an  existing  damage  model  and  failure  criterion  available  in  LS- 
DYNA  (ex.  the  Johnson-Cook  model  with  damage,  the  GISSMO  (Generalized  Incremental  Stress 
State  dependant  damage  MOdel)  in  the  *mat  add  erosion  card,  etc.).  Or  it  could  be  use  a 
particular  method  (see  the  ones  presented  in  section  3.2)  or  develop  a  new  one.  In  either  case,  a 
study  on  the  size  of  the  mesh  will  be  needed  to  predict  damage/failure  in  the  compartment. 
Experimental  data  should  also  be  used  to  validate  the  model. 

Finally,  one  alternative  to  predict  damage  on  the  compartment  would  be  to  look  at  the  stresses 
developed  in  the  structures  instead  of  eroding  elements  in  the  structure.  This  method  would  be 
simpler  and  would  not  require  the  development  and  validation  of  a  damage/failure  model. 
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4  General  Conclusion  and  Recommandations 


This  document  presented  an  investigation  on  the  feasibility  of  using  LS-DYNA  to  perform  the 
finite  element  structural  analysis  to  study  the  structural  damage  of  a  typical  compartment-type 
vessel.  The  work  was  divided  in  two  main  parts:  various  studies  were  conducted  to  generate  and 
propagate  a  free  air  blast  and  a  discussion  on  damage  and  failure  models  was  initiated. 

In  the  first  part  of  the  work,  four  studies  were  presented  to  generate  and  propagate  a  free  air  blast. 
In  each  study,  a  large  domain  that  could  contain  a  typical  compartment  was  meshed. 

The  first  and  second  studies,  demonstrated  that  an  ALE  model  with  a  rectangular  mesh  with 
element  size  of  200  mm,  100  mm,  75  mm  and  50  mm  could  not  predict  accurately  the  maximum 
incident  pressure  and  impulse  at  distances  of  1.5  m,  2.0  m  and  3.0  m  from  the  center  of 
detonation.  In  order  to  capture  accurate  shock  pressure  and  impulse,  a  rectangular  mesh  with 
element  size  in  the  order  of  1  to  1 0  mm  would  probably  be  required. 

The  third  study  presented  the  *Ioad_blast_enhanced  (LBE)  approach  combined  to  a  multi¬ 
material  arbitrary  lagrangian  eulerian  (ALE)  model  in  LS-DYNA.  The  study  showed  that  we  can 
use  the  LBE-ALE  technique  with  element  size  of  100  mm  and  get  representative  maximum 
incident  pressure  and  impulse  on  the  first  layer  of  elements  of  the  ALE  domain.  However,  only  a 
direct  loading  has  to  be  considered  (everywhere  where  there  are  no  ALE  elements,  reflections 
from  walls,  floor  or  ceiling  are  not  considered).  Also,  for  this  work,  since  CONWEP  impulse  was 
only  valid  for  target  distances  beyond  7.57  charge  radii,  threat  close  to  the  walls,  ceiling  or  floor 
cannot  be  simulate  using  the  LBE-ALE  approach. 

In  the  fourth  study,  a  spherical  mesh  was  used  to  propagate  a  free  air  blast  detonated  in  the  center 
of  the  compartment.  The  simulation  demonstrated  that  a  spherical  mesh  with  3,222,400  elements 
provided  maximum  incident  pressures  and  impulse  closer  to  the  CONWEP  predictions  compared 
to  a  rectangular  mesh  with  3,517,344  elements.  Another  advantage  of  the  spherical  mesh 
compared  to  the  rectangular  mesh  was  that  there  is  no  directional  dependence  of  the  mesh. 
However,  the  loading  resulting  from  multiple  reflections  would  not  be  accurate.  Also,  this  model 
cannot  be  used  to  simulate  a  detonation  everywhere  in  the  box  as  it  would  require  increasing  the 
sphere  radius  to  still  enclose  the  entire  compartment  when  moving  the  spherical  mesh  and 
consequently,  it  would  increase  the  number  of  elements. 

Other  methods  could  have  been  tested,  however,  due  to  the  short  time  given  to  realise  this  work, 
it  was  not  possible  to  investigate  all  of  them.  Considering  the  limitations  of  the  different 
approaches  to  simulate  an  internal  blast  it  is  suggested  to  continue  the  development  and  validation 
of  the  raytracer  approach. 

In  the  second  part  of  the  work,  the  most  important  aspects  to  consider  when  modeling 
damage/failure  of  elements  in  a  FE  (finite  element)  model  were  highlighted.  It  was  suggested  to 
develop  a  method  to  predict  the  damage  and  failure  in  large  scale  model. 

This  method  could  be  to  use  an  existing  damage  model  and  failure  criterion  available  in  LS- 
DYNA  (ex.  the  Johnson-Cook  model  with  damage,  the  GISSMO  (Generalized  Incremental  Stress 
State  dependant  damage  Model)  in  th Q^mat  add  erosion  card,  etc.  Or  it  could  be  use  a  particular 
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method  (see  the  ones  presented  in  section  3.2)  or  develop  a  new  one.  In  either  case,  a  study  on  the 
size  of  the  mesh  will  be  needed  to  predict  damage/failure  in  the  compartment.  Experimental  data 
should  also  be  used  to  validate  the  model. 

Finally,  the  last  alternative  suggested  to  predict  damage  on  the  compartment  was  to  look  at  the 
stresses  developed  in  the  structures  instead  of  eroding  elements  in  the  structure.  This  method 
would  be  simpler  and  would  not  require  the  development  and  validation  of  a  damage/failure 
model. 
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Annex  A  Range  of  Applicability  of  CONWEP 


Schwer  (2010-A)  presented  a  comparison  between  the  pressure  and  impulse  of  a  TNT  free  air 
burst  predicted  by  CONWEP  and  given  by  Swizdak  (1975)  for  4.5  kg  and  454.5  kg  charges.  The 
comparisons  between  them  showed  that  the  maximum  incident  pressures  agreed  quite  well; 
however,  CONWEP  overpredicted  the  impulses  at  close  range. 

In  his  paper,  Schwer  (2010-A)  explained  how  Kingery  and  Bulmash  (1984)  parameterized  a 
collection  of  experimental  data  using  TNT  Equivalence  of  different  explosives  and  cube  root 
scaling  of  range,  impulse  and  time  with  explosive  weight  combined  with  the  Friedlander  equation 
into  a  form  that  can  be  used  by  CONWEP.  Kingery  and  Bulmash  (1984)  provided  the  range  of 
applicability  (in  charge  radii)  for  the  positive  phase  duration  for  a  spherical  free  air:  2.7  to  750 
charge  radii.  Therefore,  usually  the  minimum  number  of  applicability  used  is  3  times  the  charge 
radii  when  using  CONWEP. 

Schwer  (2010-A)  also  gave  the  range  of  applicability  in  terms  of  scaled  distance  (Z)  by  extracting 
data  from  Swizdak  report.  The  scaled  distance  given  was  Zi=0.404  to  28.18  kg/m13  for  4.5  kg 
charge  and  Z2=0.404  to  26.18  kg/m13  for  454.5  kg  charge.  It  is  possible  to  transform  these  scaled 
distances  in  terms  of  charge  radii.  The  scaled  distance,  comes  from  the  cube  root  scaling  method 
and  is  given  by: 


Z  = 


R 

W 173 


(1) 


where  R  is  the  range  from  the  spherical  charge  (m),  W  is  the  mass  of  the  charge  (kg).  For  4.5  kg, 
using  equation  (1),  the  Ri  (range  from  the  spherical  charge)  obtained  is  0.667  m  and  for  454.5  kg, 
the  R2  obtained  is  3.106  m.  The  spherical  charge  radius  r  is  given  by  equation  (2): 


'  TNTch  arg  e 


3  W 


,1/3 


\jtp 


TNT  J 


(2) 


Where  p  is  the  density  of  TNT  and  equals  1570  kg/m3  [Kingery  and  Bulmash  (1985)]  and  W  is 
the  mass  of  TNT  (kg).  For  4.5  kg,  ri=0.0881  m  and  for  454.5  kg,  r2=0.4104  m. 

Finally,  the  charge  radii  is  given  by  R^p  =  7.57  and  R2/r2  =  7.57. 

Using  the  scaled  distance,  the  charge  radii  was  calculated  to  be  7.57  for  a  100  kg  TNT  charge. 
Therefore,  since  according  to  Schwer  (2010-A),  the  maximum  incident  pressures  agreed  quite 
well  for  the  range  of  distances  given;  a  charge  radii  of  3  can  be  used  for  comparison  with 
CONWEP’ s  predictions  but  since  CONWEP  oveipredicted  the  impulses  at  the  closest  ranges,  a 
charge  radii  of  7.57  (corresponding  to  1.852  m  for  a  100  kg  TNT  charge)  should  be  used  in  this 
work  to  compare  the  impulse  data  with  CONWEP’s  predictions. 
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Annex  B  Incident  pressure  versus  time  examples 


This  annex  presents  a  comparison  between  the  incident  pressure  histories  obtained  for  a 
rectangular  mesh  at  distances  of  1.5  m,  2.0  m  and  3.0  m  from  the  center  of  detonation  for  three 
rectangular  meshes  with  element  sizes  of  100  mm,  75  mm  and  50  mm. 

Figure  C-l:  Incident  pressure  for  a  rectangular  mesh  with  element  size  of  100  mm  at  1.5  m,  2.0  m 

and  3. 0  m. 
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Figure  C-2:  Incident  pressure  for  a  rectangular  mesh  with  element  size  of  75  mm  at  1.5  m,  2.0  m 

and  3. 0  m. 


Time  (ms) 


Figure  C-3:  Incident  pressure  for  a  rectangular  mesh  with  element  size  of  50  mm  at  1.5  m,  2.0  m 

and  3. 0  m. 
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Annex  C  GVAM/IN BLAST  State  of  Work 


This  annex  summarizes  the  state  of  work  of  GV AM/INBLAST  from  July  25th,  2012. 

1.  Stephane  modified  the  INBLAST93  code  to  read  GVAM  INBLAST  files.  INBLAST93 
includes  a  raytracer  for  shock  modeling  as  well  as  the  former  QS  pressure  propagation 
algorithms.  INBLAST93  had  no  response  algorithms:  (QS)  failure  pressures  and  minimum 
loading  durations  are  user- input.  These  are  normally  predicted  by  the  GVMFM1  pre¬ 
processor  based  on  web,  frame  and  stiffener  spacing  and  resulting  panel  eigenfrequencies. 
Note  that  a  lognormal  damage  algorithm  is  used  to  account  for  incontrollable  variations  in 
panel  toughness  and  loading;  the  output  rupture  pressure  is  below  the  nominal  one 
(conservative),  although  a  true  randomisation  would  be  better  (several  INBLAST  runs 
would  be  averaged) 

2.  RIPTASP  was  implemented  in  this  for  modeling  rigid-plastic  cross-stiffened  panel 
response,  as  a  replacement  to  the  GVAM  elastic  response  algorithms  not  trusted  anymore. 
Note  that  RIPTAB  beam  algorithms  are  automatically  used  by  RIPTASP  if  there  is  a  single 
row  of  sub-panel  cells.  A  limitation  is  that  the  code  assumes  a  spatially  uniform  loading. 

3.  Another  possibility  would  have  been  to  use  a  Dynamic  Load  Factor  in  conjunction  with  a 
(probably  clamped  panel)  rupture  criterion  as  in  Gass  et  al.  (1988)  (EXBLAST),  but  this 
did  not  work  well  in  EXBLAST.  Note  that  RIPTAB  was  already  implemented  as  an 
alternative  in  EXBLAST. 

4.  (QS)  Rupture  for  RIPTASP  is  to  be  based  on  the  maximal  strain  in  plastic  parts  of  the  sub¬ 
beams.  The  user  will  decide  what  level  constitutes  a  rupture  (2%  is  the  onset  of,  25%  is 
quasi  certain  rupture).  Note  that  GVAM  INBLAST  also  lets  the  user  decide  the  hole  size  in 
case  of  rupture. 

5.  The  above  works  in  quasi-static  (QS)  cases.  For  shocks,  since  RIPTAB  assumes  a  spatially 
uniform  loading,  it  was  suggested  using  the  above  approach,  but  after  defining  an 
equivalent  uniform  pressure  obtained  by  equating  the  associated  moments  M  on  a  clamped 
edge  associated  to  a  spatially  uniform  and  a  distributed  load,  respectively.  This  implies  a 
clumsy  recalculation  at  each  time  step  and  all  nodes.  Another  approach  would  be  to  use  an 
impulse  criterion  as  P.  27  of  the  above  report  or  Chapter  3  of  TM5-1300. 

The  above  formula  is: 


I  C  =  SD  *  h  *  A*  B  * 


(i  +  aQ(i-2aQ 
e{i-n) 


(3) 


Where  I_C=critical  impulse  =  !(P(x,y,t)  dt  dx  dy)  of  positive  pressure  history  over  whole 
plate  (N*s) 
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SD  =  Dynamic  yield  strength  of  material  (110%  of  static  one  for  steel  following  old  TM5- 
855  9-2)  (Pa) 
h  =  plate  thickness  (m) 

A,  B  =  in-plate  dimensions  (m) 
rho  =  plate  density  (kg/m3) 

N  =  Poisson  ratio 
E  =  Young’s  modulus  (Pa) 

This  is  simply  based  on  the  usual  shock  propagation  formula  by  putting  S  =  SD.  S  is  given 
by  the  actual  stress  (Pa): 


S  =  rho  *C*V 


(4) 


V  is  the  plate  material  velocity  (here  averaged  in  space  coordinates)  (m/s) 
C  is  the  in-plate  shock  velocity  (m/s)  given  by: 


rho  * 


(l  +  /V)(l  — 2/V) 
E{\  -  N) 


(5) 


6.  The  above  approach  applies  obviously  to  single  shocks  only.  If  a  reflected  shock  impinges 
the  panel  before  rebound  (if  any)  this  will  not  work.  The  report  suggests  a  (n  obscure) 
method  for  dealing  with  this  situation.  Note  that  even  if  the  panel  has  time  to  rebound 
before  arrival  of  a  reflected  shock,  it  may  nevertheless  be  deformed  permanently,  which 
may  affect  the  yield  pressure,  thus  affecting  the  effect  of  the  second  shock,  another  subtlety 
not  accounted  for  by  the  above  method. 

7.  Ideally,  the  raytracer  should  give  a  history  load  to  the  nodes  in  LS-DYNA.  This  could 
require  some  simplifications. 
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Young’s  Modulus  (Pa) 
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